Polariton crystallization in driven arrays of lossy nonlinear resonators 
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We investigate the steady states of a lossy array of nonlinear optical resonators that are driven by 
lasers and interact via mutual photon tunneling. For weak nonlinearities, we find two-mode squeez- 
ing of polaritons in modes whose quasi-momenta match the relative phases of the laser drives. For 
strong nonlinearities the spatial polariton density-density correlations indicate that the polaritons 
crystallize and are predominantly found at a specific distance from each other despite being injected 
by a coherent light source and damped by the environment. 
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Introduction - Interacting quantum many-body sys- 
tems [l| give rise to a number of fascinating phenomena 
such as quantum phase transitions, quantum magnetism 
or charge fractionalization. Of particular interest is the 
strongly correlated regime, where collective phenomena 
are most pronounced. In most cases however, strongly 
correlated quantum many-body systems are studied in 
scenarios of thermodynamic equilibrium, that permit a 
description with statistical techniques, and substantial 
understanding of, e.g. equilibrium quantum phase tran- 
sitions Q, has been achieved. On the other hand, a lot 
less is known about non-equilibrium regimes where the 
balance between loading and loss mechanisms leads to 
the emergence of stationary states. Here we investigate 
collective phenomena in non-equilibrium steady states of 
lossy arrays of coupled nonlinear optical resonators that 
are coherently driven by lasers. 

Strongly interacting polaritons 0, H[ and photons in 
coupled arrays of micro-cavities @ and optical fibers [H, 
0] have recently been shown to be suitable candidates for 
realizing a strongly correlated many-body regime with 
current technology [1,0]. So far, possibilities to observe 
equilibrium phenomena, such as a Mott insulator 0, [H[ or 
a Tonks-Girardeau gas jg,0| , have mostly been addressed. 
These regimes have however been realized previously, e.g. 
with ultra cold atoms [l|. In contrast, we here predict 
a phenomenon for polaritons, for which no analogue in 
other implementations is known so far. 

Experiments to generate quantum states with photons, 
either in cavity QED [i| or with optical fibers Q, typi- 
cally work in non-equilibrium situations and it is there- 
fore much more natural and feasible to consider driven 
dissipative scenarios. First steps in this direction have 
been undertaken with studies of an optical Josephson ef- 
fect (Io| . the dynamical evolution for nonlinearities ini- 
tially prepared in a non-equilibrium state [IH , an analy- 
sis of the spectroscopical properties of driven dissipative 
nonlinearities and entanglement studies fl3j ] . 

In this work, we consider arrays of cavities that are 
driven by lasers of constant intensity and dissipate pho- 
tons into their environment. Photons can tunnel between 
neighboring cavities and interact with suitable emitters 



in each cavity in such a way that they form polaritons 
and experience an optical nonlincarity. In this scenario, 
the interplay of laser drive and photon loss leads to the 
emergence of steady states, for which we derive the par- 
ticle statistics and characteristic correlations. 

We find two main results. In the regime, where the 
Rabi frequencies of the driving lasers are much stronger 
than the nonlinearities, only one Bloch mode with quasi- 
momentum k is driven by the lasers, and we find two- 
mode squeezing for modes with quasi-momenta p and p, 
such that p + p = k. Since this sqeezing emerges for 
weak nonlinearities, an experimental observation would 
not require a strong coupling regime for the employed 
cavities. In the complementary regime, where the nonlin- 
earities are much stronger than the laser drives, we find 
spatial anti-correlations of the polariton densities indi- 
cating that polaritons crystallize and are predominantly 
found within a specific distance from one another. As 
it requires strong nonlinearities, the crystallization can 
only be generated in cavities that operate in the strong 
coupling regime Q. We stress that this polariton crys- 
tallization appears despite the fact that coherent lasers 
continuously drive the cavities and damping permanently 
dissipates photons. The emergence of crystallization in 
the dissipative scenario with coherent drive is the most 
significant result of this work and has no analogue in 
other realizations. 

Model - Since bare photons do not interact, photon- 
photon interactions or optical nonlinearities only emerge 
when light interacts with optical emitters. Depending 
on the strength of the photon-emitter coupling, the ele- 
mentary excitations of the system are either photons, for 
weak coupling, or polaritons, superpositions of photons 
and emitter excitations, for strong coupling. In the fol- 
lowing we will use the term "polaritons" for both regimes. 
For both regimes, their Hamiltonian can be taken to read, 
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in a frame that rotates at the frequency uil of the driving 
lasers (we set h = 1). We assume periodic boundary con- 
ditions, the index j labels the resonators and aj (a,-) cre- 
ates (annihilates) a polariton in resonator j. Polaritons 
interact with strength U in each resonator and tunnel be- 
tween neighboring resonators at rate J. A = uj po i — wt is 
the detuning between polariton and laser frequency and 
Q,j are the Rabi frequencies of the driving lasers. We as- 
sume that all lasers have the same amplitude, but may 
have different phases, f2j = ile~ 4< ^ . Only relative phases 
of the lasers matter and we can choose f2 > 0. For now, 
we choose = %j (j = 1, 2, . . . , N) and TV to be a mul- 
tiple of 4 for reasons that will become clear in the sequel. 
Other values of cj>j will be considered below. 

The Hamiltonian ([1} can be implemented in several 
ways One suitable approach Q makes use of dark 
state polaritons in 4-level atoms, where a dispersively op- 
erated two polariton process gives rise to the interaction 
term Y a ] a j( a ] a j — 1) m each resonator. The dynamics 
of the system, including polariton losses from the cavities 
at a rate 7 is given by the master equation 
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FIG. 1: The steady state in the strong driving regime. Left: 
n as given by eq ([3]) as a function of Q/7 and U/j. Right: 
m/n as given by eqs. ([5]) and |3| as a function of fi/7 and 



and for the steady state, the density of photons in the 
background field, n = |/3| 2 , is determined by the equa- 
tion 4U 2 n 3 + 7 2 ?i = Q 2 , which has 



n = (3 1 / 3 * 2 / 3 - 3 2 / 3 7 2 )/(6C/^ 1/3 ) 
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The Hamiltonian H of eq 
ten in terms of Bloch modes, B k 
where k = 
to read H 
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arbitrary integer z and uj^ = A — 2J cos A:. The damping 
terms transform to Ylki^^kP^l ~ bIb^p — pBtB^). 

For our specific choice of TV and 4>j — lasers that 
drive each cavity resonantly (A = 0), constructively in- 
terfere in driving the mode Bjl of the same frequency 
ujiL = A = 0. Lasers that are in phase, <pj = 4>q, would 
destructively interfere for this mode, B * , thus motivating 
our choice of N and VI j . We note that the lasers generate 
a polariton flow in the cavity array, that can roughly be 
estimated as ~ Jsin0, where <j> = iln^fij+i/fij) is phase 
difference between the driving lasers of adjacent cavities, 
and becomes maximal for <j) = tt/2. We now analyze the 
steady states of eq. (O, for which p = 0. 

Strong driving regimes - The k = ir/2 mode is, in 
contrast to all other modes driven by the lasers and po- 
laritons from this mode can only scatter into other modes 
via the nonlinearities U. For regimes, where D, U, 
one thus expects that the state of the polariton field in 
the cavity array can be well approximated by a coherent 
state in the mode k = tt/2 plus small perturbations. We 
therefore split the mode operators, B % = /3k + b k , into 
coherent parts, represented by a complex number and 
quantum fluctuations bk, where j3k = VN /35k, Neglect- 
ing all quantum fluctuations, bk, the background field /3 



obeys the equation of motion, j5 
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with X = 9Un 2 + V3VV + 27c7 2 £l 4 as the only real and 
positive solution. Furthermore, a stability analysis (l4| 
shows, that this solution is always stable, which guar- 
antees the existence of a unique steady state. The left 
plot in figure [1] shows n as a function of E//7 an d fi/7. 
n is maximal for U = and O 3> 7. Expanding the 
Hamiltonian to second order in bk and b\ we obtain, 
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where k = -Mt-k — k and terms that are linear in b?_ have 
|fc| _ 2 

been neglected as they cancel in the corresponding mas- 
ter equation by virtue of eq. ([3]). The Hamiltonian (0| 
is known to lead to two-mode squeezing for the pairs of 
modes (k, k) (l7j . The corresponding master equation 
(|2|) with H as in eq. d3| is quadratic in the operators 
b\ and bk- Its steady state is therefore a Gaussian state, 
that is completely determined by the first and second 
order moments of b* k and bk, which are zero except for, 
{b\bk) = m and (bkbk) = g, where 
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We can now check the validity of our approximation by 
verifying that X)fc(^I^fc) <C TVn m <C n. The result- 
ing phase diagram is shown in the right plot of figure [TJ 
where we plot m/n as a function of C//7 and fi/7. In 
the regime with f2 > 7, there is on average more than 
one polariton in each cavity, n > 1, and for increasing 
nonlinearities U , the state differs significantly from a co- 
herent state in mode k = tt/2. For SI < 7, on the other 
hand, the polariton density is small, n < 1. Since the 
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nonlincarities only affect states with more than one po- 
lariton, they become ineffective in this regime and the 
state remains coherent for higher values of U. 

To obtain a more detailed picture of the steady state 
we study its particle statistics, which can be calcu- 
lated via its characteristic function 171. For the first 



order coherence between modes k and p we obtain 
Nn + Sk, p m, whereas the density- 
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density correlations between modes k and p, gm (k,p) = 
(B{BlB p B k )/((BiB k )(BtB p )), read, g^\k,p) « 1 + 

<5fc/§ (Sk, P + 4+ P ,±7r^2-) - S k,% s p,%m to leading or- 
der in l/(Nn). Here we have taken into account that 
\g\ 2 ,m 2 <C n 2 . Due 
the driven mode k = 

ffm^f,!) « 1 - jr The nonlinearity U always 



to the weak nonlinearity £/, 
? shows slight anti-bunching 



scatters two polaritons from mode k = ^ synchronously 
into modes k and k, c.f. eq. (U]). This leads to two- 
mode squeezing where polariton pairs as given by g in 

i i 2 

cq. (O are created. The term Sk+ P ,±n-^r describes cor- 

relations that originate from these pairs. Since = 
4 + 7 2 /(4n 2 £/ 2 ), these are most pronounced for strong 
damping, U/j -C 1 and O/7 -C 1, where however the 
intensity becomes increasingly weak, n C 1, sec fig. [1] 
Importantly the observation of these pairing correlations 
does not require a strong coupling regime for the em- 
ployed cavities. Nonetheless, the effect is not classical 
and disappears for U — > 0, where g — > 0, c.f. eq. (0. 
For the correlations between different resonators, j 

and I, g[ l \j,l) = (a]ai) / ^ (a]aj)(a\ai) and g r (2) (j',0 = 

{a\a\aiaj)/({a\aj){a\ai)), we find for \g\ 2 ,m 2 
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# « e^U-'J, g^(jj) » 1 - 2f and hence, 

as expected for a nonlinearity, a spatial anti-bunching of 
gP(jJ) «l-2* 0. 

It is important to note, that all above results are in- 
dependent of the tunneling rate J, a feature that only 
appears for our specific choice, 4>j = ^j, for the phases of 
the lasers. For laser phases where \A—Jcos(<f>j—<f>j-i)\ > 
^7, bistabilities of the steady state appear. A detailed 
discussion of these cases will be presented elsewhere. 

Weak driving regimes - We now analyze the regime, 
where the Rabi frequencies of the driving lasers fij are 
weaker than the nonlinearities U. For these parameters 
we have to resort to numerical calculations. We represent 
the density matrix of the polaritons as a Matrix Prod- 
uct Operator and employ a TEBD algorithm [l6| that 
integrates the master equation ([2]) in time until a steady 
state is reached. We use a 2nd order Trotter decompo- 
sition for the Lindblad super-operator with time steps 
St = [//100 (or St = [7/50) and keep relative errors due 
to matrix truncation below 10 -8 (or 10~ 6 ) at each time 
step allowing for matrix dimension up to 300 x 300. We 
consider an array of 16 cavities and, since O -C U, trun- 
cate the Hilbert space to allow for up to 2 polaritons in 




FIG. 2: Density correlations of the steady state for 16 cavities 
in the weak driving regime with A = 0, U/j = 10, fi/7 = 
2 and J/7 = 2. Top left: g?\j,l). Top right: g£ ] (k,p). 
Bottom left: g? ] (8, 1) and 3^(8, 1). Bottom right: gr (2) (8,Z) 
and p^ 2 q(8, I), zoomed in on 8 < j < 16. 



each cavity. 

In a first example we choose A = 0, U/j = 10, 
0/7 = 2 and J/7 = 1. Figure [2] shows density correla- 



tions, gr 2 \j, I) (top left) and gm'(k,p) (top right), for the 
steady state. As a consequence of the strong nonline arity 
there is a pronounced anti-bunching, g^ 2 \j,j) -C 1 [17| . 
In marked contrast to the strong driving regimes, polari- 

(2) 

tons in the same mode are strongly paired, g m (k, k) > 
1, but polaritons in different modes are strongly anti- 

(21 

correlated or anti-paired, g m (k,p) < 1 for k ^ p. That 
is, if a polariton is found in mode k, the probability to 
find a second polariton in a mode p 7^ k is lower than 

(2) 

for independent particles. g m (k,p) is smallest for k ^ p, 
where one quasi-momentum equals 7r/2. 

Most interestingly, g^ 2 \j,j + 1) is larger than unity 
whereas g^{j,j + 2) and g^ 2 \j,j + 3) etc. are signifi- 
cantly below unity, see bottom row of figure [2l Polari- 
ton densities in neighboring cavities are thus correlated 
whereas they are anti-correlated for larger separations. 
That is, if a polariton is found in one cavity, the proba- 
bility to find a second polariton is for separation 1 higher 
and for larger separations lower than for independent par- 
ticles. This behavior indicates that the polaritons are 
crystallized and predominantly occur at distances of one 
cavity-cavity separation from each other. More specifi- 
cally, the polaritons form dimers that are extended across 
two neighboring resonators and move along the array due 
to the flow created by the relative phases of the lasers. 

Since our system is one-dimensional, has significantly 
less than one polariton per cavity and polaritons in the 
same cavity strongly interact, one might be tempted to 
compare it to a Tonks-Girardeau gas [l[. For the lat- 
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FIG. 3: gi'(8,j) and n(8) as functions of J/7. Left: ^'(8,8) 
and n(8) as functions of J/7. Right g^(8,j) for 8 < j < 16 
and J/7 = 0, 0.1, 0.5, 1 and 2. A = 0, C//7 = 10 and fi/7 = 2. 
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FIG. 4: Left: gl 2) (8,j) for A = 0, J/7 = 1, fi/7 = 2 and 
C//7 = 5. Right: ^ 2) (8, j) for A = 0, t7/ 7 = 10, J/7 = 1 and 
n/7 = 2 for = 0,?r/4 and tt/2. 

ter, g^ 2 ) is the same as for free fermions and shows os- 
cillating anti-correlations, which are known as Friedel 
oscillations [HI and appear as a consequence of the 
Pauli exclusion principle. The density anti-correlations 

we find are quantitatively different from Friedel oscilla- 

(2) 

tions of free Fermi or Tonks-Girardeau gases, S^gU' = 
l- r^y ) 2 , where n is the number of polaritons per 

cavity. This is shown in the bottom row of fig. where 

(2) (2) 
we plot gi (8, j) and 5tq(8,j) with h = n(8). In con- 
trast to gi , g^Q only shows anti-correlations, g^c — 
Whereas the amplitude of the anti-correlations is compa- 
rable, they do not oscillate. 

Next we consider the dependence of the densities, h, 
and density correlations, g^ 2 \j,l), on the polariton tun- 
neling J in more detail. We have computed g^ 2 '(8,j) 
and n(8) for J/7 = 0,0.1,0.5,1 and 2, where the other 
parameters are A = 0, U/ r ) — 10 and O/7 = 2. The 

left plot of figure [3] shows gi- (8,8) and n(8) as func- 

(2) 

tions of J/7, whereas the right plot shows g? (8,j) for 
8 < j < 16 and J/7 = 0, 0.1, 0.5, 1 and 2. The crystalliza- 
tion signatures appear for nonzero tunneling J only and 
density anti-correlations become increasingly pronounced 
and long ranged as J is increased. 

To confirm the experimental robustness of our findings 

(2) 

we have computed <?r (8, j) for smaller nonlinearities, U, 
and various phase differences, <p = i ln(J7j+i/f2j). Figure 
H shows g? ) {8 1 i) for A = 0, J/7 = 1 and fi/7 = 2. In 
the left plot we choose <fi = ir/2 and set C//7 = 5. In the 
right plot we choose U/j = 10 and consider <f> = 0,7r/4 
and 7r/2. Anti-correlations only appear for <fi =/= 0. 



(2) 

Whereas the strong anti-bunching, gr ^ 1, 

is expected for a strong nonlinearity [13], the anti- 

(2) 

correlations, g} < 1 for \j — l\ > 2, are more sur- 

prising. They emerge due to and interplay between the 
nonlinearities and the polariton flow generated by the 
relative phase differences of the lasers, which in turn is 
scattered at the nonlinearities. 

In all our examples we find h(j) < 0.5 and g^ 2 \j,j) <SC 
1, which confirms that truncating the local Hilbert space 
to states of at most 2 polaritons is indeed a good ap- 
proximation. The validity of our calculations is also sub- 
stantiated by their excellent agreement with the exact 
analytical solution for the J = limit (bij ]. 

Experimental realization and measurements - The 
crystallization of polaritons we predicted here can be 
observed with resonators of high single emitter cooper- 
ativity, such as microtoroids, circuit cavities, photonic 
band gap cavities, micropillar Bragg stacks or Fabry- 
Perot microcavities on a silicon chip |8|. A straight for- 
ward method to measure the correlations we derived is 
to detect the light emitted from the structure. Detection 
of near-field photons with detectors of sufficiently fast re- 
sponse time gives access to correlations between cavities, 
ffr (J, I) and g[ 2 \j,l) whereas the far-field carries infor- 
mation on correlations between the modes, gm(k,p) and 
(2) 

<7m (k,p). Furthermore, the polaritons arc superpositions 
of photons and emitter excitations and their statistics 
and correlations can be inferred from measurements on 
the emitters. In some implementations, the polaritons 
can even be perfectly transferred onto the emitters prior 
to the measurement 0]. Even though the variations of 

(2) (2) 

#r (j 7^ are only in the range 0.95 < g^ (j 7^ I) < 

(2) 

1.05, they can reliably be measured since c/r (j 7^ I) is a 
ratio of density correlations which are both affected by 
detector inefficiencies in the same way, leaving their ratio 
unaffected. 

The author thanks M. Kiffner and A. Recati for discus- 
sions. This work is part of the Emmy Noether project HA 
5593/1-1 funded by Deutsche Forschungsgemeinschaft 
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